Średnia modelowana jako proces AR: "AR" (np. AR(1), AR(2), …)
Ocena wydajności modelu
Testowanie istotności parametrów
Test t-Studenta
Dla każdego parametru \(\theta_i\) modelu GARCH przeprowadzamy test istotności, aby sprawdzić, czy jego wartość jest statystycznie różna od zera. W tym celu używa się testu t-Studenta, który sprawdza hipotezę zerową \[
H_0: \theta_i = 0,
\] przeciwko \[
H_1: \theta_i \neq 0.
\] Statystyką testową jest \[
t_i = \frac{\hat{\theta_i}}{SE(\hat{\theta_i})},
\] gdzie \(\hat{\theta_i}\) to oszacowana wartość parametru, a \(SE(\hat{\theta_i})\) to jego błąd standardowy. Asymptotycznie \(t_i\) ma rozkład normalny, więc wykonujemy sprawdzenie \[
|t_i| > z_{1-\alpha/2},
\] gdzie \(z_{1-\alpha/2}\) to kwantyl rozkładu normalnego standardowego dla poziomu istotności \(\alpha\).
Alternatywnie, możemy patrzeć na poziom krytyczny \(p\)-value, który jest obliczany jako \[
p = 2 \cdot (1 - \Phi(|t_i|)),
\] gdzie \(\Phi\) to funkcja dystrybuanty rozkładu normalnego. Jeśli \(p < \alpha\), to odrzucamy hipotezę zerową i uznajemy, że parametr jest istotny statystycznie.
W praktyce, zazwyczaj, gdy \(|t_i| > 2\), to możemy uznać, że parametr jest istotny na poziomie \(\alpha = 0.05\).
Testowanie istotności parametrów
Istotność parametrów dla modelu GARCH(1,1)
import numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom arch import arch_model
Wczytujemy dane i obliczamy logarytmiczne stopy zwrotu z akcji Google. Wykres logarytmicznych stóp zwrotu przedstawia zmienność cen akcji w czasie.
aapl = pd.read_csv('aapl.us.txt', index_col=0)aapl['return'] = aapl['Close'].pct_change()aapl = aapl['return']aapl.index = pd.to_datetime(aapl.index, format='%Y-%m-%d')aapl.dropna(inplace=True)aapl.plot(title='Returns of Apple Stock', figsize=(14, 4))
Dopasowujemy model.
model = arch_model(aapl, p =1, q =1, mean ='constant', vol ='GARCH', dist ='t')model_fit = model.fit(disp='off', cov_type ='robust')model_fit.summary()
Wiemy, że \(r_t|\mathcal{F}_{t-1} \sim N(\mu, \sigma_t^2)\), czyli mozemy napisać funkcję gęstości prawdopodobieństwa: \[
f(r_t|\mathcal{F}_{t-1};\theta) = \frac{1}{\sqrt{2\pi\sigma_t^2}} \exp\left(-\frac{(r_t - \mu)^2}{2\sigma_t^2}\right),
\] gdzie \(\theta = (\mu, \omega, \alpha, \beta)\) to wektor parametrów modelu. Wtedy funkcja wiarygodności to iloczyn funkcji gęstości prawdopodobieństwa dla wszystkich obserwacji: \[
L(\theta|\textbf{r}) = \prod_{t=1}^{T} \frac{1}{\sqrt{2\pi\sigma_t^2}} \exp\left(-\frac{(r_t - \mu)^2}{2\sigma_t^2}\right).
\] Logarytmujemy funkcję wiarygodności, aby uzyskać funkcję log-wiarygodności: \[
\ln L(\theta|\textbf{r}) = -\frac{T}{2} \log(2\pi) - \frac{1}{2} \sum_{t=1}^{T} \ln(\sigma_t^2) - \frac{1}{2} \sum_{t=1}^{T} \frac{(r_t - \mu)^2}{\sigma_t^2}.
\] Wartości \(\sigma_t^2\) są wyznaczane rekurencyjnie z modelu GARCH(1,1), czyli np. dla \(t=1\) mamy \[
\sigma_1^2 = \omega + \alpha \cdot \epsilon_{0}^2 + \beta \cdot \sigma_{0}^2,
\] dla \(t=2\) mamy \[
\sigma_2^2 = \omega + \alpha \cdot \epsilon_{1}^2 + \beta \cdot \sigma_{1}^2 = \omega + \alpha \cdot (r_1 - \mu)^2 + \beta \cdot \sigma_{1}^2.
\] Jako warunki początkowe wybiera się \(\sigma_0^2 = \frac{\omega}{1-\alpha-\beta}\) albo \(\sigma_0^2 = \text{wariancja próbkowa }r_t\) oraz \(\epsilon_0 = 0\).
Czyli maksymalizujemy wektor \[
S(\theta|\textbf{r}) = \begin{pmatrix}
\frac{\partial \ln L}{\partial \mu} \\ \frac{\partial \ln L}{\partial \omega} \\ \frac{\partial \ln L}{\partial \alpha} \\ \frac{\partial \ln L}{\partial \beta}
\end{pmatrix}.
\] Następnie, aby znaleźć wartości parametrów, które maksymalizują funkcję wiarygodności, musimy użyć numeryki, ponieważ nie ma prostego analitycznego rozwiązania.
Odwracana macierz to macierz informacji Fishera (ujemna macierz Hessego).
Jest to asymptotyczny estymator macierzy kowariancji parametrów modelu.
Zakłada, że model i założony rozkład residuów są poprawne.
Robust Covariance Matrix (Bollerslev-Wooldridge)
W przypadku modelu GARCH(1,1) z rozkładem normalnym, możemy użyć estymacji błędów standardowych Bollerslev-Wooldridge’a, które są bardziej odporne na heteroskedastyczność i autokorelację w resztach modelu. Ta macierz kowariancji ma postać \[
\text{Cov}_{\text{BW}}(\hat{\theta}) = (-\textbf{H})^{-1} \cdot \textbf{J} \cdot (-\textbf{H})^{-1},
\] gdzie \(\textbf{H}\) to macierz Hessego zdefiniowana jak wyżej, a \(\textbf{J} = \sum_{t=1}^T \frac{\partial \ell}{\partial \theta} \cdot \frac{\partial \ell}{\partial \theta}^T\) i \(\ell_t(\theta) = \frac{1}{2}\ln(2\pi)-\frac{1}{2}\ln\sigma_t^2-\frac{1}{2}\frac{(r_t-\mu)^2}{\sigma_t^2}\) to logarytm wiarygodności w punkcie \(t\).
Jeśli residua są leptokutyczne albo mają warunkową heteroskedastyczność, gradienty będą uwzględniały dodatkową zmienność.
Ta metoda jest przydatna jeśli residua nie są gaussowskie.
Obliczenie błędów standardowych
Szuakane błędy standardowe estymowanych parametrów to pierwiastki kolejnych elementów na przekątnej macierzy kowariancji, czyli \[
\text{SE}(\hat{\theta_i}) = \sqrt{\text{Cov}(\hat{\theta})_{ii}}.
\]
Walidacja założeń modelu
Wykres residuów
Szukamy wzorców w resztach modelu, które mogą sugerować, że model nie jest odpowiedni. Wykresy residuów powinny być losowe i nie wykazywać żadnych wyraźnych wzorców.
Zgodnie z założeniami reszty powinny być z rozkładu normalnego, widać jednak często występujące ekstremalne wartości.
from statsmodels.graphics.tsaplots import plot_acffig, ax = plt.subplots(1, 1, figsize=(14, 4))plot_acf(model_fit.std_resid, lags=40, ax=ax, alpha=0.05)ax.set_title('ACF of Standardized Residuals')plt.show()
Walidacja założeń modelu
Test Ljung-Boxa
Test Ljung-Boxa sprawdza, czy reszty modelu są niezależne. Hipotezy testowe są następujące: \[
H_0: \text{reszty są niezależne} \\
H_1: \text{reszty są skorelowane}.
\] Robimy wykres \(p\)-value dla różnych opóźnień \(h\).
from statsmodels.stats.diagnostic import acorr_ljungboxljung_box = acorr_ljungbox(model_fit.std_resid, return_df=True)ljung_box['lb_pvalue'].plot(title='Ljung-Box Test p-values', figsize=(14, 4), style='.')plt.axhline(y=0.05, color='r', linestyle='--')plt.xlabel('Lag')plt.ylabel('p-value')plt.show()
Miary dopasowania modelu
Akaike Information Criterion (AIC)
\[
AIC = -2 \cdot \log(L) + 2 \cdot k
\] gdzie \(L\) to funkcja wiarygodności, a \(k\) to liczba parametrów w modelu. Im mniejsza wartość AIC, tym lepsze dopasowanie modelu.
Preferuje modele o mniejszej liczbie parametrów, ale nie karze ich zbytnio za złożoność.
Jeśli zależy nam na lepszych zdolnościach predykcyjnych, to lepiej użyć AIC.
AIC dla naszego modelu wynosi -38655.03.
Bayesian Information Criterion (BIC)
\[
BIC = -2 \cdot \log(L) + k \cdot \log(n)
\] gdzie \(L\) to funkcja wiarygodności, \(k\) to liczba parametrów w modelu, a \(n\) to liczba obserwacji. Podobnie jak AIC, im mniejsza wartość BIC, tym lepsze dopasowanie modelu.
Preferuje modele o mniejszej liczbie parametrów, ale karze je bardziej za złożoność.
Jeśli zależy nam na lepszym dopasowaniu modelu do danych, to lepiej użyć BIC.
Wykonujemy prognozę o jeden krok do przodu, a następnie porównujemy prognozowane wartości z rzeczywistymi wartościami. Wykonujemy prognozę na podstawie modelu GARCH(1,1) i obliczamy MAE i RMSE.
QUASI-MAXIMUM LIKELIHOOD ESTIMATION AND INFERENCE IN DYNAMIC MODELS WITH TIME-VARYING COVARIANCES https://public.econ.duke.edu/~boller/Published_Papers/ectrev_92.pdf